!\Description:
!  This subroutine returns the converged approximations to eigenvalues
!  of A*z = lambda*B*z and (optionally):
!  There is negligible additional cost to obtain eigenvectors.  An orthonormal
!  basis is always computed.  There is an additional storage cost of n*nev
!  if both are requested (in this case a separate array Z must be supplied).
!
!  The approximate eigenvalues and eigenvectors of  A*z = lambda*B*z
!  are derived from approximate eigenvalues and eigenvectors of
!  of the linear operator OP prescribed by the MODE selection in the
!  call to ZNAUPD.  ZNAUPD must be called before this routine is called.
!  These approximate eigenvalues and vectors are commonly called Ritz
!  values and Ritz vectors respectively.  They are referred to as such
!  in the comments that follow.   The computed orthonormal basis for the
!  invariant subspace corresponding to these Ritz values is referred to as a
!  Schur basis.
!\Arguments:
!  SELECT  Logical array of dimension NCV.  (INPUT)
!          If HOWMNY = 'S', SELECT specifies the Ritz vectors to be
!          computed. To select the  Ritz vector corresponding to a
!          Ritz value D(j), SELECT(j) must be set to .TRUE..
!          If HOWMNY = 'A' or 'P', SELECT need not be initialized
!          but it is used as internal workspace.
!
!  D       complex(8) array of dimension NEV+1.  (OUTPUT)
!          On exit, D contains the  Ritz  approximations
!          to the eigenvalues lambda for A*z = lambda*B*z.
!
!  Z       complex(8) N by NEV array    (OUTPUT)
!          On exit, if RVEC = .TRUE. and HOWMNY = 'A', then the columns of
!          Z represents approximate eigenvectors (Ritz vectors) corresponding
!          to the NCONV=IPARAM(5) Ritz values for eigensystem
!          A*z = lambda*B*z.
!
!          If RVEC = .FALSE. or HOWMNY = 'P', then Z is NOT REFERENCED.
!
!          NOTE: If if RVEC = .TRUE. and a Schur basis is not required,
!          the array Z may be set equal to first NEV+1 columns of the Arnoldi
!          basis array V computed by ZNAUPD.  In this case the Arnoldi basis
!          will be destroyed and overwritten with the eigenvector basis.
!
!  LDZ     Integer.  (INPUT)
!          The leading dimension of the array Z.  If Ritz vectors are
!          desired, then  LDZ .ge.  max( 1, N ) is required.
!          In any case,  LDZ .ge. 1 is required.
!
!  SIGMA   complex(8)  (INPUT)
!          If IPARAM(7) = 3 then SIGMA represents the shift.
!          Not referenced if IPARAM(7) = 1 or 2.
!
!  WORKEV  complex(8) work array of dimension 2*NCV.  (WORKSPACE)
!
!  **** The remaining arguments MUST be the same as for the   ****
!  **** call to ZNAUPD that was just completed.               ****
!
!  NOTE: The remaining arguments
!
!           BMAT, N, WHICH, NEV, TOL, RESID, NCV, V, LDV, IPARAM, IPNTR,
!           WORKD, WORKL, LWORKL, RWORK, INFO
!
!         must be passed directly to ZNEUPD following the last call
!         to ZNAUPD.  These arguments MUST NOT BE MODIFIED between
!         the the last call to ZNAUPD and the call to ZNEUPD.
!
!  Three of these parameters (V, WORKL and INFO) are also output parameters:
!
!  V       complex(8) N by NCV array.  (INPUT/OUTPUT)
!
!          Upon INPUT: the NCV columns of V contain the Arnoldi basis
!                      vectors for OP as constructed by ZNAUPD .
!
!          Upon OUTPUT: If RVEC = .TRUE. the first NCONV=IPARAM(5) columns
!                       contain approximate Schur vectors that span the
!                       desired invariant subspace.
!
!          NOTE: If the array Z has been set equal to first NEV+1 columns
!          of the array V and RVEC=.TRUE. and HOWMNY= 'A', then the
!          Arnoldi basis held by V has been overwritten by the desired
!          Ritz vectors.  If a separate array Z has been passed then
!          the first NCONV=IPARAM(5) columns of V will contain approximate
!          Schur vectors that span the desired invariant subspace.
!
!  WORKL   Double precision work array of length LWORKL.  (OUTPUT/WORKSPACE)
!          WORKL(1:ncv*ncv+2*ncv) contains information obtained in
!          znaupd.  They are not changed by zneupd.
!          WORKL(ncv*ncv+2*ncv+1:3*ncv*ncv+4*ncv) holds the
!          untransformed Ritz values, the untransformed error estimates of
!          the Ritz values, the upper triangular matrix for H, and the
!          associated matrix representation of the invariant subspace for H.
!
!          Note: IPNTR(9:13) contains the pointer into WORKL for addresses
!          of the above information computed by zneupd.
!          -------------------------------------------------------------
!          IPNTR(9):  pointer to the NCV RITZ values of the
!                     original system.
!          IPNTR(10): Not used
!          IPNTR(11): pointer to the NCV corresponding error estimates.
!          IPNTR(12): pointer to the NCV by NCV upper triangular
!                     Schur matrix for H.
!          IPNTR(13): pointer to the NCV by NCV matrix of eigenvectors
!                     of the upper Hessenberg matrix H. Only referenced by
!                     zneupd if RVEC = .TRUE. See Remark 2 below.
!          -------------------------------------------------------------
!
!  INFO    Integer.  (OUTPUT)
!          Error flag on output.
!          =  0: Normal exit.
!
!          =  1: The Schur form computed by LAPACK routine csheqr
!                could not be reordered by LAPACK routine ztrsen.
!                Re-enter subroutine zneupd with IPARAM(5)=NCV and
!                increase the size of the array D to have
!                dimension at least dimension NCV and allocate at least NCV
!                columns for Z. NOTE: Not necessary if Z and V share
!                the same space. Please notify the authors if this error
!                occurs.
!
!          = -1: N must be positive.
!          = -2: NEV must be positive.
!          = -3: NCV-NEV >= 1 and less than or equal to N.
!          = -5: WHICH must be one of 'LM', 'SM', 'LR', 'SR', 'LI', 'SI'
!          = -6: BMAT must be one of 'I' or 'G'.
!          = -7: Length of private work WORKL array is not sufficient.
!          = -8: Error return from LAPACK eigenvalue calculation.
!                This should never happened.
!          = -9: Error return from calculation of eigenvectors.
!                Informational error from LAPACK routine ztrevc.
!          = -10: IPARAM(7) must be 1,2,3
!          = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatible.
!          = -12: HOWMNY = 'S' not yet implemented
!          = -13: HOWMNY must be one of 'A' or 'P' if RVEC = .true.
!          = -14: ZNAUPD did not find any eigenvalues to sufficient
!                 accuracy.
!          = -15: ZNEUPD got a different count of the number of converged
!                 Ritz values than ZNAUPD got.  This indicates the user
!                 probably made an error in passing data from ZNAUPD to
!                 ZNEUPD or that the data was modified before entering ZNEUPD
!-----------------------------------------------------------------------
subroutine zneupd(nev, d, z, n, sigma, which, tol, resid, ncv, iparam, ipntr, workd, workl, lworkl, info)
    implicit none
    character  which*2
    integer    info, n, lworkl, ncv, nev
    complex(8) sigma
    real(8) tol
    ! Array Arguments
    integer iparam(11), ipntr(14)
    logical select(ncv)
    real(8) rwork(ncv)
    complex(8) d(nev), resid(n), z(n, nev), workd(3 * n), workl(lworkl)
    ! Local Scalars
    complex(8), parameter :: one = (1.0D0, 0.0D0), zero = (0.0D0, 0.0D0)
    complex(8) rnorm, temp, vl(1), workev(2 * ncv)
    real(8) conds, sep, rtemp, eps23
    integer bounds, ierr, ihbds, iheig, nconv, invsub, iuptri, j, mode, irz, ibd, outncv, numcnv, k
    logical reord
    ! External Functions
    real(8), external :: dznrm2, dlamch, dlapy2
    complex(8), external :: zdotc

    ! Executable Statements
    mode = iparam(7)
    nconv = iparam(5)
    info = 0

    ! Get machine dependent constant.
    eps23 = dlamch('Epsilon-Machine')
    eps23 = eps23**(2.0D+0 / 3.0D+0)
    ! Quick return Check for incompatible input
    ierr = 0
    if (nconv <= 0) then
        ierr = -14
    else if (n <= 0) then
        ierr = -1
    else if (nev <= 0) then
        ierr = -2
    else if (ncv <= nev .or.  ncv > n) then
        ierr = -3
    else if (which /= 'LM' .and. which /= 'SM' .and. which /= 'LR' .and. which /= 'SR' .and. &
             which /= 'LI' .and. which /= 'SI') then
        ierr = -5
    else if (lworkl < 3 * ncv**2 + 4 * ncv) then
        ierr = -7
    end if

    if (ierr /= 0) then
        info = ierr
        return
    end if
    !--------------------------------------------------------%
    ! Pointer into WORKL for address of H, RITZ, WORKEV, Q   |
    ! etc... and the remaining workspace.                    |
    ! Also update pointer to be used on output.              |
    ! Memory is laid out as follows:                         |
    ! workl(1:ncv*ncv) := generated Hessenberg matrix        |
    ! workl(ncv*ncv+1:ncv*ncv+ncv) := ritz values            |
    ! workl(ncv*ncv+ncv+1:ncv*ncv+2*ncv) := error bounds     |
    !--------------------------------------------------------%

    ! The following is used and set by ZNEUPD. 
    ! workl(ncv*ncv+2*ncv+1:ncv*ncv+3*ncv) := The untransformed Ritz values.
    ! workl(ncv*ncv+3*ncv+1:ncv*ncv+4*ncv) := The untransformed error bounds of the Ritz values
    ! workl(ncv*ncv+4*ncv+1:2*ncv*ncv+4*ncv) := Holds the upper triangular matrix for H.
    ! workl(2*ncv*ncv+4*ncv+1:3*ncv*ncv+4*ncv) := Holds the associated matrix representation of
    !                                      the invariant subspace for H.
    ! GRAND total of NCV * ( 3 * NCV + 4 ) locations.
    bounds = ipntr(8)
    iheig = bounds + ncv
    ihbds = iheig + ncv
    iuptri = ihbds + ncv
    invsub = iuptri + ncv * ncv

    ! irz points to the Ritz values computed by _neigh before exiting _naup2.
    ! ibd points to the Ritz estimates computed by _neigh before exiting _naup2.
    irz = ipntr(14) + ncv * ncv
    ibd = irz + ncv
    ! RNORM is B-norm of the RESID(1:N).
    rnorm = workl(3)
    workl(3) = 0.0D0

    reord = .FALSE.
    ! Use the temporary bounds array to store indices. These will be used to mark the select array later
    select(:) = .FALSE.
    do j = 1, ncv
        workl(bounds + j - 1) = j
    end do
    ! Select the wanted Ritz values.
    ! Sort the Ritz values so that the wanted ones appear at the tailing
    ! NEV positions of workl(irr) and workl(iri).  Move the corresponding
    ! error estimates in workl(ibd) accordingly.
    call zsortc(which, ncv, workl(irz), workl(bounds))
    ! Record indices of the converged wanted Ritz values
    ! Mark the select array for possible reordering
    numcnv = 0
    do j = 1, ncv
        rtemp = max(eps23, dlapy2(dble(workl(irz + ncv - j)), dimag(workl(irz + ncv - j))))
        k = workl(bounds + ncv - j)
        if (numcnv < nconv .and. dlapy2(dble(workl(ibd + k - 1)), dimag(workl(ibd + k - 1))) <= tol * rtemp) then
            select(k) = .true.
            numcnv = numcnv + 1
            if (k > nev) reord = .true.
        end if
    end do

    !-----------------------------------------------------------%
    ! Check the count (numcnv) of converged Ritz values with    |
    ! the number (nconv) reported by dnaupd.  If these two      |
    ! are different then there has probably been an error       |
    ! caused by incorrect passing of the dnaupd data.           |
    !-----------------------------------------------------------%
    if (numcnv /= nconv) then
        info = -15
        return
    end if

    !-------------------------------------------------------%
    ! Call LAPACK routine zlahqr to compute the Schur form  |
    ! of the upper Hessenberg matrix returned by ZNAUPD.    |
    ! Make a copy of the upper Hessenberg matrix.           |
    ! Initialize the Schur vector matrix Q to the identity. |
    !-------------------------------------------------------%
    call zcopy(ncv * ncv, workl, 1, workl(iuptri), 1)
    call zlaset('All', ncv, ncv, zero, one, workl(invsub), ncv)
    call zlahqr(.true., .true., ncv, 1, ncv, workl(iuptri), ncv, workl(iheig), 1, ncv, workl(invsub), ncv, ierr)
    call zcopy(ncv, workl(invsub + ncv - 1), ncv, workl(ihbds), 1)

    if (ierr /= 0) then
        info = -8
        return
    end if

    if (reord) then
        ! Reorder the computed upper triangular matrix.
        call ztrsen('None', 'V', select, ncv, workl(iuptri), ncv, workl(invsub), ncv, workl(iheig), &
                nconv, conds, sep, workev, ncv, ierr)

        if (ierr == 1) then
            info = 1
            return
        end if
    end if

    ! Copy the last row of the Schur basis matrix to workl(ihbds).  This vector will be used
    ! to compute the Ritz estimates of converged Ritz values.
    call zcopy(ncv, workl(invsub + ncv - 1), ncv, workl(ihbds), 1)
    !--------------------------------------------%
    ! Place the computed eigenvalues of H into D |
    ! if a spectral transformation was not used. |
    !--------------------------------------------%
    if (3 /= mode) then
        call zcopy(nconv, workl(iheig), 1, d, 1)
    end if
    !----------------------------------------------------------%
    ! Compute the QR factorization of the matrix representing  |
    ! the wanted invariant subspace located in the first NCONV |
    ! columns of workl(invsub,ldq).                            |
    !----------------------------------------------------------%
    call zgeqr2(ncv, nconv, workl(invsub), ncv, workev, workev(ncv + 1), ierr)

    !--------------------------------------------------------%
    ! * Postmultiply V by Q using zunm2r.                    |
    ! * Copy the first NCONV columns of VQ into Z.           |
    ! * Postmultiply Z by R.                                 |
    ! The N by NCONV matrix Z is now a matrix representation |
    ! of the approximate invariant subspace associated with  |
    ! the Ritz values in workl(iheig). The first NCONV       |
    ! columns of V are now approximate Schur vectors         |
    ! associated with the upper triangular matrix of order   |
    ! NCONV in workl(iuptri).                                |
    !--------------------------------------------------------%
    call zunm2r('R', 'N', n, ncv, nconv, workl(invsub), ncv, workev, z, n, workd(n + 1), ierr)

    do j = 1, nconv
        !---------------------------------------------------%
        ! Perform both a column and row scaling if the      |
        ! diagonal element of workl(invsub,ldq) is negative |
        ! I'm lazy and don't take advantage of the upper    |
        ! triangular form of workl(iuptri,ldq).             |
        ! Note that since Q is orthogonal, R is a diagonal  |
        ! matrix consisting of plus or minus ones.          |
        !---------------------------------------------------%
        if (dble(workl(invsub + (j - 1) * ncv + j - 1)) < 0.0D0) then
            call zscal(nconv, -one, workl(iuptri + j - 1), ncv)
            call zscal(nconv, -one, workl(iuptri + (j - 1) * ncv), 1)
        end if
    end do

    ! Compute the NCONV wanted eigenvectors of T located in workl(iuptri,ldq).
    do j = 1, ncv
        if (j <= nconv) then
            select(j) = .true.
        else
            select(j) = .false.
        end if
    end do

    call ztrevc('R', 'S', select, ncv, workl(iuptri), ncv, vl, 1, workl(invsub), ncv, ncv, outncv, workev, rwork, ierr)

    if (ierr /= 0) then
        info = -9
        return
    end if

    ! Scale the returning eigenvectors so that their Euclidean norms are all one. LAPACK subroutine
    ! ztrevc returns each eigenvector normalized so that the element of largest magnitude has magnitude 1.
    do j = 1, nconv
        rtemp = dznrm2(ncv, workl(invsub + (j - 1) * ncv), 1)
        rtemp = 1.0D0 / rtemp
        call zdscal(ncv, rtemp, workl(invsub + (j - 1) * ncv), 1)

        ! Ritz estimates can be obtained by taking the inner product of the last row of the
        ! Schur basis of H with eigenvectors of T.
        ! Note that the eigenvector matrix of T is upper triangular, thus the length of the
        ! inner product can be set to j.
        workev(j) = zdotc(j, workl(ihbds), 1, workl(invsub + (j - 1) * ncv), 1)
    end do

    ! Copy Ritz estimates into workl(ihbds)
    call zcopy(nconv, workev, 1, workl(ihbds), 1)
    ! The eigenvector matrix Q of T is triangular.  Form Z*Q.
    call ztrmm('R', 'U', 'No', 'Non', n, nconv, one, workl(invsub), ncv, z, n)

    ! Transform the Ritz values and possibly vectors and corresponding error bounds of OP to those of A*x = lambda*B*x.

    ! A spectral transformation was used. Determine the Ritz estimates of the Ritz values in the original system.
    call zscal(ncv, rnorm, workl(ihbds), 1)
    if (3 == mode) then
        do j = 1, ncv
            temp = workl(iheig + j - 1)
            workl(ihbds + j - 1) = workl(ihbds + j - 1) / temp / temp
        end do
        !-----------------------------------------------------------%
        ! *  Transform the Ritz values back to the original system. |
        !    For TYPE = 'SHIFTI' the transformation is              |
        !             lambda = 1/theta + sigma                      |
        ! NOTES:                                                    |
        ! *The Ritz vectors are not affected by the transformation. |
        !-----------------------------------------------------------%
        do j = 1, nconv
            d(j) = 1.0D0 / workl(iheig + j - 1) + sigma
        end do
        ! Eigenvector Purification step. Formally perform one of inverse subspace iteration.
        !------------------------------------------------%
        ! Purify the computed Ritz vectors by adding a   |
        ! little bit of the residual vector:             |
        !                      T                         |
        !          resid(:)*( e    s ) / theta           |
        !                      NCV                       |
        ! where H s = s theta.                           |
        !------------------------------------------------%
        do j = 1, nconv
            if (workl(iheig + j - 1) /= zero) then
                workev(j) = workl(invsub + (j - 1) * ncv + ncv - 1) / workl(iheig + j - 1)
            end if
        end do
        ! Perform a rank one update to Z and purify all the Ritz vectors together.
        call zgeru(n, nconv, one, resid, 1, workev, 1, z, n)
    end if
end
